Fast Fourier Transform Simulation Techniques for Coulomb Gases 
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An improved approach to updating the electric field in simulations of Coulomb gases using the 
local lattice technique introduced by Maggs and Rossetto Ij, is described and tested. Using the 
Fast Fourier Transform (FFT) an independent configuration of electric fields subject to Gauss' 
law constraint can be generated in a single update step. This FFT based method is shown to 
outperform previous approaches to updating the electric field in the simulation of a basic test 
problem in electrostatics of strongly correlated systems. 
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I. INTRODUCTION 

Electrostatic forces play a major role in intra- and 
inter-molecular interactions. In particular, in locally 
charged systems, strong Coulomb interactions between 
charged particles often dominate structural and dynam- 
ical properties . This situation occurs frequently in 
biological systeiiis, where both protein 0, Q and DNA 
macromolecules typically undergo H-|- dissociation or 
association, thereby generating a charged macroion with 
counterions in solution. In addition, aqueous solutions 
often contain dissolved salts, thus introducing more mo- 
bile ions into the system. The thermally averaged in- 
teraction force (which integrates to give the "potential 
of mean force" (PMF) ^, JJ) between two macroions de- 
pends strongly on electrostatic forces that are supplied by 
the mobile ions in the surrounding solution. Calculation 
of the PMF, or, equivalently, electrostatic free energy of 
interaction between two macroions, is critical to under- 
standing their equilibrium (and, ultimately, dynamical) 
properties. For example, the electrostatic PMF deter- 
mines the binding equilibrium constant for the macroion 
pair 0. Similar remarks hold for synthetic macroions, for 
example, charged colloid particles, whose self-assembly 
properties, phase diagrams, etc., are determined by pre- 
cisely the same type of considerations "9, 10]. Dynamical 
properties, such as rate constants for transitions between 
stable states and ion currents through channel pro- 
teins prj , are heavily influenced by the energy landscape 
determined from equilibrium free energy considerations. 

The difficulty of implementing such long range ( "non- 
local") forces in realistic numerical simulations is well 
known: each charged particle "feels" the effect of all 
the others, so that naive algorithms typically scale with 
the square of the system volume V , while more sophisti- 
cated algorithms (Ewald summation, Fourier techniques 
[T^ I13I ll^) with better scaling properties nevertheless 
prove far from efficient in computing the electrostatic en- 



ergy for large systems. Recently, a local formulation of 
the Coulomb gas problem (first introduced by Maggs and 
collaborators |3] has offered a possible exit from this im- 
passe. At the cost of including an unphysical transverse 
(or "curl" ) part in the electric field, the effects of which on 
the charged particle dynamics can be shown to decouple 
on average over a suitably long Monte Carlo simulation 
of the system, the Hamiltonian of the system can be re- 
cast in a completely local form, so that the computational 
cost of a system update becomes of the order of the sys- 
tem volume V (rather than , as above) . This method, 
together with several improvements designed to improve 
the mobility of the charged particles in the simulations, 
have been the subject of a number of recent publications 
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The local Coulomb gas formalism is based on the parti- 
tion function of a system consisting of N charges (mobile 
or fixed) at locations fi, thereby producing a charge 
density 
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Accordingly, the canonical partition function for the sys- 
tem at inverse temperature /3 becomes 

„ N . 
Z= W dT,VE{f) \{5(y-E- —p{T))e-^ I "'^"^"^ 
i=i f ^ 

(2) 

where e is the dielectric constant of the medium in which 
the charges move (assumed spatially uniform throughout 
this paper). The essential point is that the transverse, or 
curl, part of the electric field variable decouples from the 
charged particle dynamics via the Helmholtz decomposi- 
tion 
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+ V xA (3) 
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As only the gradient part of the electric field appears in 
the Gauss' Law constraint in jSJ, the curl part clearly 
decouples from the particle positions. However, the sim- 
ulation of the integral in Q involves charged particle 
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moves that are only sensitive to the local value of the 
electric field, provided that the local field is modified to 
maintain Gauss' Law at all times. Practically, the simula- 
tions are most readily performed on spatial lattices, with 
the particle charges associated with lattice sites, and the 
electric fields with links connecting nearest neighbors on 
the lattice. For a complete description of the algorithm, 
the reader is referred to Ref. jlj. 



II. PREVIOUS WORK 

Before discussing our FFT based approach to updating 
the electric field in a local Coulomb gas simulation, we 
briefly review some of the field update methods that have 
been used in previous applications of the local lattice 
method of Maggs et al ^ JJj . 

The implementation of the method is most simply seen 
on a square lattice with dimensionless variables. We 
therefore scale the electric field by a factor of ea^/Ane and 
the inverse temperature by a factor of 47re^/ea, where e is 
the dielectric constant of the system and a is the lattice 
spacing. 

Any field update method must conserve Gauss' law at 
every step of the simulation. The simplest methods for 
updating the electric field are based on shifting the four 
link fields of a plaquette by the same amount a. Such 
a shift maintains the Gauss' law constraint. One possi- 
bility is to use a Metropolis procedure to pick a, but it 
is more efficient to use a heat-bath method. Consider- 
ing one plaquette and labeling its links 1 through 4, the 
portion of the Hamiltonian that depends on these links 
(written using dimensionless variables) is given by 

I 4 

Hp^-J2iE^ + af (5) 

1=1 

where Ei is the electric field on the ith link in the direc- 
tion obtained by going around the plaquette clockwise. 
This can be rewritten as 

H, = lia-J2E,r (6) 

i 

where we have dropped terms with no a dependence. 
Therefore a can be generated with the canonical dis- 
tribution by setting it to x/\J~^ — '^^Ei where x is a 
Gaussian-distributed variable with zero mean and stan- 
dard deviation one and /3 is the dimensionless inverse 
temperature. Updating the plaquettes using this proce- 
dure will be referred to in this paper as the local heat 
bath method, as the change in each plaquette is affected 
only by the immediate local neighborhood of that plaque- 
tte. Of the three methods compared in this paper, the 
local heat bath is the simplest to implement but also the 
"most local," and can be expected to yield the longest 
autocorrelation times. 



Level, Alet, Rottler and Maggs [T^l recently introduced 
a cluster move to update the electric fields more effi- 
ciently based on a worm algorithm developed by Alet 
and S0rensen . The method introduces a pair of posi- 
tive and negatively charged particles on a randomly cho- 
sen site of the lattice. One of the particles then moves 
about the lattice in a biased random walk, as described 
in Ref. • This pseudo-particle moves through the lat- 
tice (modifying appropriately the electric field to preserve 
Gauss' law), eventually returning to its partner where 
the two annihilate each other. A loop of links with mod- 
ified electric fields is left behind along the path taken 
by the mobile particle. This change of the electric fields 
along this loop is then accepted or rejected depending on 
the fields entering the initial site. Empirically these loop 
changes are accepted with high probability. The charge 
of the virtual pair is randomly chosen with a uniform dis- 
tribution in the interval {—em, Sm), where e™ is a free 
parameter which is chosen to maximize efficiency. 

The worm method is clearly more "global" than the 
local heat bath method described above — typically the 
number of links modified is of order the lattice volume — 
and we may expect autocorrelation times that are shorter 
than those obtained with plaquette updates. Neverthe- 
less, field configurations that differ by a single worm up- 
date are still substantially correlated with each other. 
Next, we show that a global regeneration of the electric 
field is possible, using FFT techniques, which is compu- 
tationally more efficient than the worm method on rea- 
sonably large lattices, and produces a completely decor- 
related electric field after a single update step. 

III. FAST FOURIER DECORRELATION 
METHOD 

In this section we describe in some detail a FFT-based 
approach to updating the electric field in local Coulomb 
gas simulations that is computationally of order V log V, 
and produces a globally decorrelated transverse field at 
each update. It is similar to methods used to study 
electromagnetic splittings in lattice quantum chromo- 
dynamics in Ref. The method does not have any 

free parameters, obviating the need to optimize simu- 
lation parameters, and, as we shall see below in some 
sample computations, is computationally more efficient 
(in the sense of computational effort per autocorrelation 
time of physically interesting observables) than the other 
updating schemes described previously. 

Once again, consider a lattice Coulomb gas with charge 
density pn and electric field link variables Enf_i (where n 
denotes lattice sites and ^ = 1,2,3 spatial direction). 
The partition function takes the form, for fixed charge 
particle locations (i.e. fixed p), 

J dEn^5{A^En^ - pn)e~^ (7) 

where A (resp. A) denote left (resp. right) lattice deriva- 
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tives here and below. We wish to develop an efficient pro- 
cedure for traversing the space of electric field configura- 
tions {En^i} consistent with the Gauss' Law constraint. 
A perfect decorrelation can be achieved if the new field 
has a transverse ("curl") part completely decorrelated 
with the previous transverse field: on the other hand, 
the longitudinal ( "gradient" ) part of the field is fixed by 
the Gauss' Law constraint and must be preserved by the 
update procedure. To accomplish this we must (a) ex- 
tract the longitudinal part, and (b) generate a new, and 
completely independent transverse part of the field. 

We shall work on a LxLxL lattice (volume V — L^) 
with Fourier transform fields defined as follows 

k 

^"/^ = -^E^'''"^M(fc) (9) 

k 

Note that we distinguish lattice coordinate space fields 
from their momentum space transforms by using sub- 
scripts for the former and function notation for the latter. 
Up to constant field configurations (which are separately 
simulated as explained in Ref. [1|), an arbitrary lattice 
vector field £'„^ may be decomposed 

En^ - Eli+E^l (10) 

A^£;lj^ = Pn (11) 
E^ifj, — ^fiup^i^^np (12) 

where Anp is a transverse lattice vector field satisfying 
^p^np — (as we see from Eq. ^1 the gradient part 
of Anp is absent from E^^'). Consequently, the Fourier 
transform of A^p, defined through 

k 

must satisfy 

s;Ap{k)^0, sp = 1 - (14) 

and can accordingly be written in terms of polarization 
vectors ei , 6*2 as follows 

Ap{k) = ai{k)eip{k) + a2{k)e2p{k) (15) 
-ei ^ -62^0 (16) 
el -el = (17) 

for a discrete mode k corresponding to a complex Fourier 
component (i.e. where k^ — 2Tm^/L with not all of 
ni,n2,n3 equal to or L/2). For the real modes the 
reader may easily verify that analogous formulas hold 
involving purely real vectors. From the form of the par- 
tition function given in Eq. [7| one finds that the Fourier 
coefficients ai, 02 in Eq.^jare to be generated according 
to the Gaussian weight 

Z(ai,a2)=e-^^«(l'^^l'+l'^^l') (18) 



with A(fc) = 4 sin^ (fcp). Once random polarization 
vectors ei and 62 are generated satisfying Equations 1161 
and llTI then Eq.^|yields the complete Fourier transform 
Ap{k). An inverse FFT then yields the coordinate space 
field Anp, and a lattice curl the desired transverse part 
of the electric field via Eq. [H 

This heat bath procedure clearly produces a new trans- 
verse electric field completely decorrelated from the pre- 
ceding one. As we wish to update the total electric field, 
we must also calculate the gradient part i?", which can 
then be added to the new decorrelated transverse field. 
To do this, we simply note that from Eq. ^] it follows 
that 

4'w = -;w^ (19) 

where p{k) is the Fourier transform of the charge den- 
sity Pn, which is also to be computed by FFT. Thus this 
algorithm involves (apart from the effort required to gen- 
erate the Fourier components Ap{k) as indicated above) 4 
FFT operations to generate a globally decorrelated elec- 
tric field. As we shall see below, this computational ef- 
fort compares favorably with competitive methods, such 
as local heat-bath or worm type update algorithms. 



IV. RESULTS 

To test our FFT-based method for updating the elec- 
tric fields, and compare its efficiency with that of the local 
heat-bath and worm update methods, we simulated the 
system of charged conducting plates with ions between 
the plates discussed in Ref. fill- The basic system con- 
sists of a 50x50x50 lattice with periodic boundary con- 
ditions in all three dimensions. Positive charges are free 
to move on two fixed plates separated in the x direction 
that extend the entire extent of the lattice in the y and 
z directions, while the region between the plates con- 
tains mobile counterions ensuring overall neutrality. As 
in Ref. the lattice spacing is chosen to be 1 A, the 
dielectric constant is 80.0, and the temperature is 30QK 
so that the dimensionless inverse temperature, /3, is 87.1. 
The plates are charged with 34 positively charged univa- 
lent ions, which are constrained to stay on the plate. 34 
negatively charged divalent ions are placed between the 
plates, and are excluded from the plane of lattice sites 
closest to the plates. Every lattice site can be occupied 
by at most one ion. The ions are moved by the coupled 
particle-field heat bath method discussed in Ref. 15]. 

Each Monte Carlo step is composed of (200 x number 
of charges on the plates) attempted moves of the positive 
ions on the plates, (2000 x number of ions in solution) 
attempted moves of negative ions in solution, a global 
update of the total electric field (as described in Ref. Ql), 
and an update of the electric field using either the local 
heat bath method, the worm method, or the FFT based 
method. For the local heat bath method 3 x 50^ links are 
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FIG. 1: Autocorrelation function of the electric field in the y 
direction at the (25,25,25) lattice site from simulations using 
FFT-based update method, local heat bath update method, 
and worm update method. After equilibration the particles 
are frozen in place for these runs. 



FIG. 2: Autocorrelation function of the (2, 2, 2) component 
of the Fourier Transform of the electric field in the x direc- 
tion from simulations using FFT-based update method, local 
heat bath update method, and worm update method. After 
equilibration the particles are frozen in place for these runs. 



chosen at random to be updated. For the worm method, 
5 worms are created for each Monte Carlo sweep. The 
maximum ghost charge used in the worm method was set 
to 0.3, which was determined experimentally to be the 
value that minimized the autocorrelation times of our 
observables. For the FFT method the FFTW package 
[20j was used to perform the FFTs. 

On a 2GHz AMD processor it takes 336 seconds to 
perform 1000 updates of the electric field using the FFT 
method. It takes 432 seconds to perform 1000 Monte 
Carlo sweeps using the local heat bath method. It takes 
700 seconds to perform 1000 Monte Carlo sweeps using 
the worm method, where each Monte Carlo sweep corre- 
sponds to creating 5 worms. Thus, using these parame- 
ters, the FFT-based method is more than a factor of 2 
faster than the worm-based method. 

To compare the efficiencies of the different methods, we 
consider the autocorrelation times of various observables. 
The autocorrelation of an observable A is given by 

where Ai designates the «th measurement of A and A is 
the average value of A. We extract the autocorrelation 
time of an observable, r, by integrating the autocorrela- 
tion function out to a distance where the measurements 
are decorrelated (t = /jJ^J"" C{t) where C{tmax) ~ 0). 

Initially we consider the autocorrelation of the elec- 
tric field on a single link of the lattice. In order to ob- 
serve how the electric field decorrelates itself using solely 
the different field update methods, we fix the location 
of the ions after the warmup sweeps. We then simulate 
the electric field using the different field update methods 
and observe the electric field in the y direction on the 
(25,25,25) site. The autocorrelation function of this elec- 



tric field is shown in Fig. ^ We determine the autocor- 
relation time by integrating the autocorrelation function 
out to 40. This gives an autocorrelation time of 7.04 for 
the worm method and an autocorrelation time of 1.61 
for the local heat bath method. The FFT based method 
gives a completely decorrelated set of electric fields, so 
the electric field under observation is decorrelated after 
a single Monte Carlo step and the autcorrelation time is 
1 . The local heat bath method is able to decorrelate this 
observable rapidly because the observable is so local. 

We next consider the autocorrelation function of the 
(2, 2, 2) component of the Fourier Transform of the elec- 
tric field in the x direction. Again we fix the location of 
the ions to concentrate on the efficiency of the field up- 
date methods. The autocorrelation functions are shown 
in Fig. 121 Integrating the autocorrelation function out 
to 50 gives an autocorrelation time of 13.3 for the local 
heat bath update method runs and 5.24 for the worm up- 
date method runs. The local heat bath method has more 
difficulties with smaller momentum components of the 
Fourier-transformed electric field: the (0, 0, 1) compo- 
nent has an autocorrelation time of 123.3 with the local 
heat bath update, while the autocorrelation time with the 
worm update is only 3.36. For this observable the large 
clusters that are updated with the worm method allow 
it to decorrelate these small-momentum Fourier compo- 
nents more rapidly then the local heat bath method. The 
FFT based method decorrelates all the Fourier compo- 
nents of the electric field in a single Monte Carlo update 
step. 

We then perform the simulations with the ions al- 
lowed to move throughout the simulation and consider 
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FIG. 3: Autocorrelation function of the particle-particle en- 
ergy from simulations using FFT based update method, local 
heat bath update method, and worm update method. 



the particle-particle energy given by the expression 
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700, we obtain an autocorrelation time of 89.85 for the 
local heat bath method simulation, 78.10 for the worm 
update method, and 47.30 for the FFT based update 
method. 



V. CONCLUSIONS 

Efficient simulation of physical systems using Monte 
Carlo simulation depends on Monte Carlo update moves 
that decorrelate the system rapidly. We have studied 
the correlations of various observables in simulations of 
charged particles using the method introduced by Maggs 
et al., and were able to reduce the autocorrelation time 
of many observables by introducing a method based on 
the FFT that produces a completely independent electric 
field (subject to the Gauss' law constraint) after a single 
update step. We have shown on a physically realistic sys- 
tem that, although local heat bath methods and worm- 
based methods are able to decorrelate the electric field 
well on different length scales, the FFT based method 
performs better than either of them on all observables 
studied. Additionally, using our implementations, the 
FFT method took the least amount of computer time. 



This is a physically interesting observable that depends 
only on the locations of the ions and not directly on 
the electric field. The correlations present in the electric 
field will affect the correlations of the particle-particle en- 
ergy. All runs are composed of 5,000 Monte Carlo equi- 
libration steps followed by 400,000 measurement steps. 
The autocorrelation functions of the particle-particle en- 
ergy for the different field update methods are shown in 
Fig. O Integrating the autocorrelation function out to 
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